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SECTION  1 
INTRODUCTION 


After  the  detonation  of  a  low  altitude  nuclear  weapon  a  strong 
blast  wave  moves  out  from  the  explosion,  eventually  (for  all  the  weapons  so 
far  detonated)  becoming  weak  and  propagating  as  an  acoustic  pulse.  As  the 
pulse  propagates  into  the  decreasing  density  of  the  upper  atmosphere  the 
pulse  strengthens  (in  terms  of  the  relative  overpressure)  and  forms  a  shock 
(more  or  less  strong  according  to  the  weapon  yield)  which  propagates  through 
l*  ionosphere.  After  the  passage  of  the  strong  pulse  the  atmosphere  and 
ionosphere  are  left  in  motion,  bobbing  up  and  down  due  to  the  interaction  of 
gravity  with  the  stratified  density  of  the  atmosphere. 

At  distances  of  thousands  of  kilomters  from  a  large,  low  altitude 
explosion,  the  ionosphere  has  been  observed  to  be  set  in  motion.  This  motion 
is  universally  attributed  to  some  type  of  internal  gravity  wave.  These 
ionospheric  motions  can  be  expected  to  produce  substantial  modifications  to 
the  propagation  of  HF  radio  waves  which  reflect  off  the  ionosphere.  In 
this  report  we  describe  a  model  for  the  ionospheric  motions  which  has  been 
programmed  into  a  ray  tracing  computer  code  which  will,  in  the  future,  be 
used  to  analyze  implications  of  gravity  wave  phenomena  for  HF  military 
systems. 


For  distances  far  from  the  burst  we  do  not  know,  in  detail,  what 
mechanisms  produce  the  gravity  wave  and  cannot  do  a  first  principles  calcu¬ 
lation  of  the  gravity  wave  effects.  For  this  region  we  must  rely  on  data 
combined  with  general  theoretical  knowledge  of  the  properties  of  gravity 
waves  to  construct  the  model.  The  model  created  for  distances  from  the  burst 
greater  than  600  km  is  de'cribed  in  Chapter  2. 
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For  distances  closer  to  the  burst  we  believe  that  our  understand¬ 
ing  of  the  important  phenomena  is  more  complete  and  we  can  perform  calcu¬ 
lations  of  the  effects  to  be  expected.  In  Chapter  3  we  describe  the  pro¬ 
cedures  by  which  we  make  these  calculations  and  give  the  results.  In  this 
chapter  we  also  describe  a  model  which  agrees  reasonably  well  with  both  the 
long  distance  model  and  the  results  of  the  calculations  for  points  closer 
in.  This  is  the  model  we  have  put  into  the  ray  trace  code. 

In  Chapter  4  we  give  the  results  of  several  examples  using  the 
model.  We  first  give  only  the  ionospheric  disturbance  calculated  by  the 
model  for  an  increasingly  complicated  set  of  inputs.  Then  we  show  the 
results  of  tracing  rays  through  the  convoluted  ionospheres  given  by  the 
model.  These  calculations  are  not  intended  to  analyze  the  systems  implica¬ 
tions  of  the  gravity  waves  but  are  intended  as  a  preliminary  demonstration 
that  the  model  works  as  expected.  We  hope  to  soon  be  able  to  use  these 
tools  to  perform  studies  in  sufficient  detail  to  make  statements  about  the 
systems  imp: ic it ions. 
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SECTION  2 


THE  MODEL  FOR  RANGES  GREATER  THAN  600  KILOMETERS 


INTRODUCTION 

The  linear  response  of  the  atmosphere  to  mild  disturbances  can  be 
given  in  terms  of  acoustic  gravity  waves  (AGW) .  Although  the  disturbance 
in  the  region  of  a  nuclear  burst  is  far  from  mild,  once  the  disturbance 
propagates  to  large  distances  its  amplitude  will  diminish  and  thus  become 
describable  in  terms  of  AGW. 

The  general  procedure  for  application  of  AGW  theory  requires  that 
the  source  be  decomposed  onto  all  possible  modes,  both  discrete  and  continuous. 
The  time  progression  of  the  system  is  then  obtained  by  advancing  the  time 
in  the  Fourier  factor  for  the  individual  modes  followed  by  a  summation  over 
the  modes.  This  procedure  can  be  reduced  to  simple  terms  in  the  event  that 
a  single  AGW  mode  describes  the  disturbance  adequately.  We  will  use  such 
simplifications.  The  approximations  used  in  modeling  ionospheric  effects 
of  AGW  are  based  on  both  theoretical  and  experimental  observations. 

GENERAL  FEATURES  OF  AGW 

In  terms  of  the  problem  at  hand,  that  of  effects  on  the  distant 
ionosphere  due  to  low-altitude  bursts,  there  is  a  set  of  theoretical  fea¬ 
tures  of  AGW  which  is  useful  to  consider.  Some  of  these  are: 

1.  The  atmosphere  supports  both  freely  propagating (continuous 
spectrum)  and  ducted  (discrete  spectrum)  modes.1  Both  may 
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propagate  to  great  distances  but  the  ducted  modes  have  their 
energy  concentrated  at  altitudes  above  50  km. 

2.  Long  period,  freely  propagating  gravity  modes  dominate  the 
distant  disturbances  launched  from  low  altitude  sources.2 
These  modes  give  a  fairly  constant  period  fluctuation  at 

a  fixed  range. 

3.  The  long  period  waves  tend  to  bend  around  the  earth's  sur¬ 
face.  3 

4.  The  particle  motion  associated  with  long  period  modes  tends 
to  be  nearly  horizontal. 

5.  The  intensity  of  the  AGW  disturbances  changes  with  altitude. 
It  tends  to  grow  exponentially  with  altitude  up  to  ~300  km 
(due  to  exponential  decrease  in  density)  and  then  decreases 
above  300  km  due  to  viscous  dissipation. 

6.  The  ionospheric  plasma  driven  by  atmospheric  motion  is 
largely  confined  to  motion  along  the  magnetic  field  lines. 
There  have  been  two  large  yield  bursts  which  have  produced 
significant  ionospheric  disturbances,  the  Pacific  test 
Housatonic,4  and  the  USSR  Novaya  Zemlya  test  of  30  October 
1961.  For  these  events,  several  features  were  consistently 
seen  by  ionosonde  observations: 

a.  The  disturbance  showed  a  nearly  sinusoidal  variation 
with  one  or  two  cycles. 

b.  These  variations  showed  a  nearly  constant  period  for 
individual  stations  which  increased  with  increasing 
distances  from  the  burst. 

c.  The  onset  time  for  the  disturbances  was  approximately 
(range) /(sound  speed)  where  the  sound  speed  is  that 
in  the  region  of  the  ionosphere. 


.-1., 
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d.  The  amplitude  of  the  disturbances  did  not  decrease  with 
distance  in  any  obvious  manner. 

e.  The  initial  motion  of  the  ionosphere  (up  or  down)  seemed 
to  be  fixed  by  the  magnetic  field  orientation. 


THE  MODEL 


In  order  to  form  a  simple  model  of  distant  ionospheric  disturb¬ 
ances  we  first  note  that  for  sinusoidal  variation  in  the  velocity  implied 
by  a  single  mode,  the  continuity  condition  is  given  by 

^(inN)  *-?•  V-V/X  (1) 

so 

£n  (N/No)  ~  V/Xo)  (2) 

where  N  is  the  electron  density  and  Nq  the  unperturbed  electron  density. 
The  position  will  be 

X  -  tQ  =J$ dt~V/w  (3) 

where  X  is  a  wavelength  and  u>  the  frequency.  Therefore  one  expects  the 
position  and  £n(N/NQ)  to  vary  in  a  related  way. 

Due  to  the  nearly  horizontal  stratification  of  the  ambient  iono¬ 
sphere,  chiefly  the  vertical  (Z)  motion  will  be  detectable.  Therefore  given 
the  magnetic  field  dip  angle  0,  and  the  azimuthal  propagation  angle  <J>, 
measured  from  magnetic  north,  purely  horizontal  (X)  neutral  motion  will  lend 
to  a  plasma  velocity  of 

vP  =  V  sin(20)  cos(<j>)/2  .  (4) 
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The  compression  implied  by  V*V,  again  confining  the  plasma  to  motion  along 
the  field  line,  gives 

V*V^  =  V  sin(20)  cos(<|))/2Az  .  (5) 

Here  A_  <<  A  has  been  used.  This  is  consistent  with  the  properties  of 

La  A 

long-period  gravity  waves. 

The  period  of  distant,  freely  propagating  gravity  waves  is 

Tc~TB/tan(a)  (6) 

where  x^  is  the  Brunt  period  and  a  the  takeoff  angle  in  Figure  1.  Thus, 
for  small  angles  a, 

TC  w  tbx/z  ' 


Figure  1.  Burst  geometry. 
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Therefore  the  argument  of  the  sinusoidal  function  describing  the  variation 
can  be  taken  as 

ARG  =  (t  -  (8) 

tbX  \  c  / 

and  this  will  satisfy  Equation  7  at  the  time  of  onset  of  the  disturbance, 

t  =  X/c  ,  (9) 

o 

where  c  is  the  speed  of  propagation. 

Although  not  addressed  specifically  to  this  point,  studies  have 
indicated  that  the  disturbance  is  nearly  independent  of  the  burst  altitude^ 
provided  the  burst  is  below  ~30  km.  These  studies  also  suggest  that  the 
effect  goes  like  the  square  root  of  the  yield5’6 and  we  will  adopt  this  con¬ 
dition.  This  dependence  is  not  well  established. 

The  intensity  of  the  disturbance  varies  with  the  altitude  of  obser¬ 
vation,  increasing  up  to  ionospheric  levels,  and  decreasing  above  the  F-region 
due  to  dissipation  from  viscous  effects.  Based  loosely  on  the  work  of  Francis1 
we  take  the  dependence  of  the  fluid  velocity  on  observation  altitude  to  be 


0  otherwise.  (10) 

Here  Z  is  the  altitude  of  peak  effect  (~300  km)  and  the  Z.  is  the 
P  L 

lowest  affected  altitude  C — 100  km). 


1-  ; 
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If  we  assume  that  the  variation  of  Equation  10  is  the  largest 
contributor  to  V*V  we  can  use  Equation  1  to  write 


d  2(z-y 

Dt  *n  N  =  - *4 

<vy 


(ii) 


which  integrates  to  give 


2 (2-2  ) 

N  =  N  Exp  — — E__ 

°  (Z,-2  )2 

L  p 


AZ 


1  - 


cz^r 


(12) 


Experiment  shows  that  gravity  waves  slowly  decay  in  strength  as 
they  propagate  away  from  the  source.  The  nuclear  data  does  not  show  this 
decay  clearly,  due  chiefly  to  the  fact  that  different  inodes  (frequencies) 
with  different  amplitudes  are  observed  at  different  distances  from  the 
bursts.  We  have  therefore  relied  on  typical  results  from  ducted  modes  to 
estimate  a  spatial  decay  factor. 


exp  ( -4  X/Cg) 


(13) 


where  Cr.  is  the  earth's  circumference. 

E 

The  model  thus  assembled  gives  values  of  AN  and  AZ  as 


Q  AZ 

AN  =  exp  — — - 

AZ  =  -  C2  •  Qj  -Q2  (14) 
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where  AZ  is  the  change  in  altitude  of  plasma  initially  at  Z  and  AN  is 
the  factor  by  which  the  electron  density  initially  at  Z  -  AZ  is  multiplied 
to  obtain  the  new  electron  density  at  Z. 


The  factors  in  Equation  14  are 

Qj  =  CY^) 1/2  cos(<J>)  sin(20)  sin  ~  (t-  |)  , 

B 

^2  = 


and 


^3  - 

2 

(yv 

Our  fit  to  the  data 

requires 

Z  = 
P 

300  km  , 

xg  =  540  sec, 

ZL  = 

100  km  , 

c  =  .65  km/sec 

CE 

4*4^  km 

CN  =  .2  (MT)*17 

and  C7  = 

20  km/MT1/2 

In  applying  this  model,  a  set  of  restrictions  are  required.  First, 


mi 


exp  (-  4X/Cg) 


(15) 


(16) 


600  km  <  X 


(17) 


since  the  region  nearly  over  the  burst  is  not  well  described.  Furthermore, 
given  that  only  two  oscillations  of  the  function  provide  the  complete  re¬ 
sponse,  there  will  be  an  effect  for 
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(18) 


x  x  2Xtb 


Also,  the  model  will  apply  only  for 


100  km  <  Z  <  500  km 


and  outside  this  region  Equation  10  provides  that  there  will  be  no  effect. 
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SECTION  3 


THE  MODEL  FOR  RANGES  LESS  THAN  600  KILOMETERS 


The  model,  described  in  detail  in  Chapter  2,  for  ranges  greater 
than  600  km  must  be  extended  to  include  the  region  nearer  the  burst.  The 
model  for  large  ranges  was  based  upon  general  properties  of  acoustic  gravity 
waves  and  upon  what  data  were  available.  We  do  not  know  in  detail  what 
physical  processes  are  involved  in  generating  and  propagating  the  long 
period  ionospheric  signals  which  have  been  observed  thousands  of  kilometers 
from  the  burst  point  so  we  are  in  no  position  to  attempt  a  first  principle 
calculation  of  the  effects  and  must  rely  heavily  on  the  data.  For  the 
region  nearer  the  burst  we  believe  we  know  what  phenomena  will  be  impor¬ 
tant  and  we  can  perform  calculations  of  the  expected  effects.  We  have 
made  these  calculations  for  a  set  of  yields  and  locations.  Below  we  give 
the  results  of  the  calculations  and  present  a  model  which  fits  these  re¬ 
sults  well  enough  and  matches  smoothly  to  the  previously  described  model 
for  large  ranges. 

The  processes  we  wish  to  calculate  and  model  are  briefly  described 
as  follows.  The  blast  wave  leaves  the  point  of  the  explosion,  assumes  the 
shape  of  an  N-wave,  that  is,  a  pulse  which  goes  disconti nuously  to  some 
positive  value,  decreases  linearly  to  a  negative  value  of  equal  magnitude, 
then  goes  discontinuously  to  zero  (see  Appendix  A)  and  weakens  sufficiently 

to  allow  for  propagation  in  the  acoustic  limit.  As  the  weak  acoustic 
pulse  propagates  into  the  rapidly  decreasing  density  of  the  higher  atmos¬ 
phere  it  restrengthens,  forming  a  shock  at  high  altitudes  which  can  in 
some  cases  be  quite  strong.  After  the  strong  N-wave  portion  of  the  upwardly 
propagating  pulse  passes  a  given  point  there  will  remain  a  residual  oscilla- 
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tion  due  to  the  buoyant  interaction  of  the  stratified  atmosphere  with  the 
earth's  gravitational  field.  In  this  report  we  will  not  attempt  to  model 

the  N-wave  itself  (we  hope  to  do  so  in  the  future)  but  will  concentrate  on 
the  residual  oscillations  since  it  is  these  oscillations  which  form  the 
"gravity  wave"  in  the  region  close  to  the  burst  point. 

In  a  previous  report  the  authors  have  described  a  procedure  for 
calculating  the  properties  of  the  N-wave  as  it  propagates  upwardly  through 
the  atmosphere.  This  report  was  originally  classified  but  has  recently  been 
declassified  and  we  have  appended  a  portion  of  the  report,  which  describes 
the  phenomenology  of  the  N-wave  pulse  in  some  detail  and  specifies  our 
procedures  for  calculating  the  properties  of  the  pulse,  to  the  present 
report  as  Appendix  A.  The  computer  code  which  calculates  the  properties  of 
the  N-wave  by  the  procedures  of  Appendix  A  is  called  MODEL. 

The  propagation  of  the  pulse  to  high  altitudes  is  an  essentially 
nonlinear  process  (the  pulse  stretches  to  become  many  times  its  original 
length).  The  residual  oscillations  which  remain  after  the  passage  of  the 
pulse  will  often  be  accurately  described  by  the  linearized  equations  of 
hydrodynamics.  We  assume  this  will  always  be  the  case.  The  linearized 
equations  of  hydrodynamics  admit  a  solution  by  means  of  a  Green's  Function 
which  has  been  given  in  closed  form.  We  can  use  the  Green's  Function  to 
represent  the  upward  propagation  of  an  N-wave  pulse  and  the  associated 
residual  oscillations  by  a  contour  integral  around  a  closed  contour  in 
complex  frequency  space.  Since  the  propagation  is  controlled  by  linear 
equations  the  length  of  the  N-wave  will  not  change  as  it  propagates.  As 
with  any  linear  system,  the  strength  of  the  pulse  is  controlled  by  an 
overall  multiplicative  factor.  As  a  part  of  a  previous  research  effort 
we  developed  a  computer  program  to  calculate  the  necessary  contour  inte¬ 
gral.  The  report  of  the  research  is  unpublished,  but  the  section  deriv¬ 
ing  the  equations  and  describing  the  program,  called  CINT,  which  performs 
the  integration  is  appended  to  the  present  report  as  Appendix  B. 


The  plan  for  the  present  calculation  is  as  follows.  We  choose 
parameters  for  a  burst  and  select  a  point  in  the  ionosphere  at  which  to 
calculate  the  gravity  wave.  We  then  run  MODEL  to  determine  the  proper 
pulse  length  and  amplitude  for  the  N-wave  at  the  point  in  question.  We 
use  the  results  from  MODEL  to  set  the  pulse  length  and  amplitude  factor 
for  CINT  which  will  then  calculate  the  parameters  of  the  gravity  wave  at 
the  point  we  wish  to  study. 

As  presently  configured  our  computer  codes  calculate  only  the 
pressure.  The  current  task  requires  that  we  have  knowledge  of  velocities 
and  displacements.  In  principle  we  could  modify  the  codes  to  calculate 
the  velocities  but  we  have  not  yet  done  so.  In  both  the  acoustic  limit 
and  the  low  frequency  limit  for  a  single  frequency  we  find  the  result  that 
v z  =  cos0  c  5p/Ypo  and  vr  =  sin  0  c  6p/ypQ  where  c  is  the  speed  of 
sound,  6p  the  pressure  perturbations,  pQ  the  ambient  pressure,  y  the 
adiabatic  constant,  and  0  the  angle,  measured  from  the  vertical, 
subtended  by  a  line  drawn  from  the  burst  point  to  the  point  for  which  we 
are  performing  the  calculation.  We  believe  that  these  relations  will  in 
general  hold  to  within  the  accuracies  required  and  we  will  use  the  codes 
to  calculate  the  pressure  fluctuations  and  then  after  the  velocities  using 
these  relations. 

The  oscillations  which  follow  the  passage  of  the  N-wave  do  not 
fall  abruptly  to  zero  after  some  time  but  die  away  more  or  less  gradually 
as  a  function  of  time.  We  will  not  try  to  describe  this  process  for  our 
model,  however,  but  will  choose  a  simple  sinusoidal  function  whose  ampli¬ 
tude  and  period  we  will  select  to  approach  the  results  from  CINT  as  best 
we  can.  To  display  for  the  reader  what  has  actually  been  done  we  show  an 
example  in  Figure  2.  The  graph  shows  the  output  from  CINT  which  resulted 
from  a  run  for  300  km  altitude  directly  above  a  three  megaton  nuclear 
explosion  at  sea  level.  The  results  from  the  code  MODEL  indicate  that  the 
length  of  the  N-wave  should  be  370  seconds  and  that  the  relative 
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Figure  2.  Sp/pQ  as  a 
a  1  RT  exp 
shown. 


i  of  time  300  km  directly  above 
The  earlier  N-wave  phase  is  not 


overpressures  at  the  shock  should  be  1.7.  The  time  of  arrival  of  the  N-wave 
is  158  seconds  so  the  residual  oscil 1 1  at  ions  which  we  call  gravity  wave,  will 
occur  after  600  seconds.  The  period  is  readily  chosen  from  the  graph  to  be 
720  seconds  but  the  choice  of  amplitude  is  not  clearcut.  We  chose  the  amplitude 
to  be  .3  and  the  duration  to  be  two  periods. 

Not  all  cases  provide  such  an  easily  determined,  constant  period 

as  that  graphed  in  Figure  2.  In  Table  1  we  have  attempted  to  organize 

the  information  on  the  period  of  the  gravity  wave  as  a  function  of  location 

and  bomb  yield.  The  entry,  L,  is  the  length  of  the  N-wave  portion  of  the 

pulse  in  seconds.  tc  is  the  natural  period  of  freely  propagating  gravity 

waves  defined  in  the  previous  chapter.  The  entries,  AT,  are  the  time 

(seconds)  between  successive  zero  crossings  during  the  gravity  wave  phase, 

that  is,  after  the  N-wave  has  passed.  We  see  that  most  often  AT  is 

reasonably  close  Cro  1/2  t  .  In  the  worst  cases  the  difference  is  a 

factor  of  two  for  the  first  AT  and  much  less  for  later  ones.  The  examples 

which  show  the  greatest  deviation  of  AT  from  1/2  tc  are  those  for 

which  L  <<  1/2  t  .  Since  the  model  for  ranges  greater  than  600  km  has  an 

oscillatory  period  of  tc  it  is  a  matter  of  considerable  convenience  to 

choose  xc  for  the  period  when  the  range  is  less  than  600  km.  The  data 

in  the  table  show  that  ic  is  probably  the  best  simple  choice  we  can 

make  so  we  fix  the  period  of  the  oscillations  at  all  ranges  to  be  T  _ 

and  the  duration  of  the  effect  at  all  ranges  to  be  2t  .  Notice  that 

c 

the  duration  is  greater  at  greater  ranges. 

We  now  attempt  to  study  the  altitude  dependence  of  the  gravity 
wave  for  ranges  less  than  600  km  from  the  burst  point.  Having  chosen 
for  the  period  of  the  oscillations  we  have  a  model  for  which  the  maximum 
displacement  will  be  A/tc  where  A  is  an  amplitude  factor  which  we  chose 
from  CINT  output  as  discussed  above.  In  Table  2  we  show  the  resulting 
displacements  for  several  altitudes  at  three  different  ranges  from  a  3  MT 
surface  burst.  In  all  cases  the  displacement  is  small  at  an  altitude  of 
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Table  1 

(Entires  are  in  seconds) 


300  km 


L=264;  i  t  =560  L=23 

2  c 


600  km 


AT=283;402;462  AT=396;663;399 


L=1 91 ;  \  tc=572 


AT=400;513;545 


300  km 


36 


600  km 


300  km 


200  km 


■sa 


AT=341 ; 359 ; 369 

AT=409;471 ;501 

AT=670  ;809 ;91 1 

L=254;  i  t  =317 

2  c 

L=337 ;  \  xc=572 

L=400;  j  t  *1003 

AT=292;290;296 

AT=406;515;543 

AT=770;964 ; 1 041 

300  km 


600  km 


300  km 


L=4b0;  g  Tc=396  L=608;  \  tc=560  1=850;  \  ic=880 


AT=328  ;361 ;  362  AT=464;478;497  AT=677;824;879 


200  km 


2;^tc=317  L=514;  |  tc=572 


1 

;  2  t 


AT=289;290;295  AT=431 ;50G;546  AT=804 ;939  ;984 
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100  km.  The  displacements  are  seen  to  peak  at  about  300  km  altitude  for 
each  of  the  three  ranges.  Recall  that  the  model  for  ranges  greater  than 

2  2 

600  km  also  had  a  peak  at  300  km  altitude,  the  function  being  [l-(z-300)  /200  ]. 
In  Figure  3  we  graph  the  curve  showing  the  altitude  dependence  for  the  model 
for  large  ranges  and  show  the  points  for  the  data  given  in  the  table  where 
we  have  normalized  the  displacements  for  each  range  to  be  one  at  300  km 
altitude.  The  points  cannot  be  said  to  all  lie  on  the  curve,  but  given 
the  crude  nature  of  our  method  for  arriving  at  the  displacements  and  the 
need  for  as  simple  a  model  as  will  suffice,  it  is  adequate  and  desirable 
to  use  the  same  altitude  dependence  for  all  ranges.  We  use  the  one  given 
in  Equation  15. 


Table  2 


\.  r 

Z  \ 

0  km 

300  km 

600  km 

500 

19 

9 

17 

400 

21 

29 

39 

Y=3MT 

300 

22 

31 

81 

200 

21 

29 

43 

The  fact  that  the  period  and  altitude  dependence  of  the  disturb¬ 
ance  can  be  taken  to  be  the  same  for  ranges  less  than  600  km  as  for  ranges 
greater  than  600  km  encourages  us  to  try  a  fairly  simple  modification  of 
Equation  15  as  a  model  for  all  ranges.  We  take  the  horizontal  neutral 
displacement  to  be 

AXN  =  '  (YMT^  ^  R  CZ^2  Sin  ^  (-20-) 

2  2  1/2 

where  R  =  (Z  +X  )  and  and  Q,  are  as  in  Equations  15  and  16.  IVe 

choose  the  neutral  displacement  in  the  vertical  direction  to  be 
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20 


(21) 


“N  ■  <V)1/2  f  %  Sin  <*  -  C>) 

In  the  limit  of  small  z/x  we  find  that  AZ^  goes  to  zero  while  AX^ 
goes  to  the  neutral  motion  on  which  the  model  for  ranges  greater  than 
600  km  was  based. 


We  have  argued  that  the  period  and  the  altitude  dependence  of 
Equations  20  and  21  are  appropriate  for  our  model;  we  must  now  check  to 
see  whether  the  amplitude  of  the  displacements  given  by  Equations  20  and 
21  agree  well  enough  with  those  which  result  from  our  calculations.  In 
Table  3  we  give  the  maximum  displacements  which  result  from  choosing  an 


Table  3. 

(Entries  are  in  km) 


r 

Y  \ 

0  km 

300  km 

600  km 

Results  from 
the  Model 

30  MT 

66 

BEK 

295 

109 

3  MT 

22 

Q 

81 

35 

0.3  MT 

6 

H 

10 

11 

amplitude  based  on  output  from  CINT.  For  each  entry  the  altitude  is  300  km, 
the  altitude  of  maximum  effect.  The  last  column  shows  results  given  by 
Equations  20  and  21  CQ2  changes  by  only  a  few  percent  in  600  km  so  the 
equations  would  predict  a  nearly  constant  amplitude).  Considering  uie 
arbitrariness  involved  in  choosing  an  amplitude  from  CINT  output  we  feel 
that  the  agreement  between  Equations  20  and  21  and  the  calculations  is 
satisfactory.  The  two  very  bad  numbers  are  the  results  at  600  km  range  for 
yields  of  3  and  30  MT.  For  these  cases  CINT  is  not  really  appropriate 
anyway  for  the  disturbances  are  so  strong  that  nonlinear  effects  will 
certainly  play  a  role.  The  nonlinear  effects  will  probably  act  to  reduce 
the  displacements.  The  displacements  given  by  Equations  20  and  21  are 
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relatively  simple.  They  agree  with  Equations  15  for  large  ranges  and 
agree  reasonably  well  with  our  calculations  for  the  region  near  the  burst. 
We  will  adopt  them  as  the  basis  for  the  gravity  wave  model. 


Equations  20  and  21  give  the  components  of  the  neutral  motion. 
We  wish  to  model  the  vertical  motion  of  the  electrons.  We  define 


AZ 


Cz (YMT^ 


1/2 


Q0  Sin 


0?  (•  ■  a) : 


COS0  sin26 


-  2sin  6 


AN  =  Exp  (Q3AZ/Q2) 


where  C  ,  Q9,  Qv  4>,  and  0,  are  as  in  Equations  15  and  16.  If  the  unper- 

Z  Z 

turbed  electron  density  is  N(Z)  we  specify  the  perturbed  density  to  be 
AN-N  (Z-AZ). 

The  fact  that  the  displacement  due  to  the  gravity  wave  scales 
like  the  square  root  of  the  yield  has  implications  with  regard  to  modeling 
multiburst  effects.  If  a  region  is  simultaneously  affected  by  gravity  waves 
from  two  bursts  whose  burst  points  are  remote  from  one  another  we  would 
expect  that  the  displacements  would  simply  add  in  a  linear  way  while  the 
compressions  would  multiply.  If  this  rule  were  to  remain  true  for  very 
close  bursts  the  gravity  wave  effects  would  scale  like  the  first  power  of 
the  yield.  We  do  not  know  the  precise  resolution  of  this  apparent  paradox 
but  we  believe  that  bursts  sufficiently  close  together  such  that  the 
resulting  blast  waves  interact  before  they  become  weak  will  produce  results 
as  a  single  burst  with  the  combined  yield  (the  square  root  scaling)  while 
bursts  further  apart  will  produce  effects  which  may  be  added  linearly. 

In  the  model  for  the  computer  programs  we  add  the  effects  from  all 
the  bursts  in  a  linear  way.  If  the  bursts  are  expected  to  be  very  close 
together  in  space  and  time  it  would  be  a  good  idea  to  combine  them  into 
a  single  burst  at  the  input  stage. 
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SECTION  4 
EXAMPLES 


In  this  chapter  we  present  the  results  of  applying  the  model  to 
specific  cases.  In  Figure  4  we  show  the  results  of  a  one  megaton  explosion 
at  sea  level.  The  contours  are  contours  of  constant  electron  density  but 
we  have  used  a  highly  artificial  ionosphere  in  which  the  electron  density 
is  set  equal  to  the  altitude.  Thus  the  plots  would  show  AZ  except  for 
the  effects  of  compression.  We  may  think  of  the  plots  as  showing  the 
effective  AZ  which  would  produce  the  same  electron  density  in  incom¬ 
pressible  flow.  Probably  the  most  important  point  to  notice  about  Fig¬ 
ures  4a  and  4b  is  the  strong  effect  of  the  magnetic  field.  In  Figure  4a 
we  show  the  propagation  of  the  gravity  wave  in  the  north-south  direction 
(magnetic  coordinates).  The  location  of  the  burst  is  at  the  center  tick 
mark  of  the  figures.  To  the  right  of  the  burst  is  magnetic  north  and  in 
that  direction  the  horizontal  velocity  moves  the  electrons  down  when  the 
vertical  velocity  moves  them  up.  The  tendency  is  to  cancel  and  the  gravity 
wave  does  not  appear  large  in  the  figure.  To  the  left  of  the  bomb,  magnetic 
south,  the  horizontal  and  vertical  components  of  the  neutral  flow  move  the 
electrons  up  and  down  in  phase  and  the  result  is  a  greatly  enhanced  wave. 

For  points  sufficiently  far  north  or  south  of  the  burst  the  vertical  vel¬ 
ocity  will  become  negligible  and  the  gravity  wave  will  exhibit  north-south 
symmetry.  We  see  clearly  in  Figure  4a  that  after  three  hours  and  at  a 
distance  at  nearly  2000  km  the  asymmetry  is  still  pronounced. 

In  Figure  4b  we  show  the  propagation  of  the  gravity  wave  in  the 
east-west  direction.  In  this  case  the  horizontal  motion  of  the  neutrals 
does  not  cause  vertical  motion  of  the  electrons  so  only  the  vertical 
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Figure  4a.  Contours  of  electron  density  is  an  artificial  ionosphere 
(see  text)  along  a  north-south  line  following  a  1  MT 
nuclear  explosion  at  sea  level.  The  vertical  tic  marks 
are  spaced  5°  ( — 556  km)  apart.  The  contour  interval  is 


Contours  of  electron  density  in  an  artificial  ionosphere 
(see  text)  along  an  east-west  line.  The  tic  marks  are 
5°  (~556  km)  apart.  The  contour  interval  is  cb  km. 


Figure  4b.  (continued) 


component  of  the  neutral  flow  can  move  the  electrons  up  and  down.  The 
effects  of  the  gravity  wave  are  therefore  confined  much  more  closely  about 
the  burst  point  and  at  a  time  of  three  hours  and  a  distance  of  1800  km  the 
effects  have  nearly  disappeared. 

In  Figure  5  we  show  contour  plots  of  a  somewhat  different  type. 

The  lines  are  contours  of  constant  electron  density  at  an  altitude  of 
250  km.  Again  we  have  used  an  ionosphere  with  the  ambient  electron  den¬ 
sity  set  equal  to  the  altitude.  The  contours  are  shown  in  increments  of 
25  km  starting  with  25  km  so  the  tenth,  or  J,  contour  is  the  ambient  value  for 
the  horizontal  plane  we  have  chosen.  For  this  figure  there  are  two  surface 
bursts,  each  10  MT.  Magnetic  north  is  to  the  right  and  down  in  the  figures. 
The  field  and  the  interaction  of  the  two  devices  with  one  another  produce 
a  complex  pattern  of  ionospheric  disturbances  even  in  this  relatively  simple 
case. 

In  Figure  6  we  show  another  view  of  the  same  two  bursts;  in  this 
case  in  a  vertical  plane  located  along  the  line  which  forms  the  X-axis 
for  Figure  5.  Along  this  line,  which  is  close  to  the  line  joining  the 
burst  points,  but  slightly  displaced,  the  ionosphere  is  seen  to  be  severely 
disturbed  for  a  period  of  hours. 

We  now  wish  to  briefly  consider  effects  on  a  realistic  communica¬ 
tions  link  from  a  portion  of  a  large  nuclear  attack.  The  HF  link,  shown  by 
the  line  in  Figure  7,  is  from  Omaha,  Nebraska  running  some  2800  km  northwest 
into  British  Columbia.  We  have  considered  100  MT  explosions  at  the  starrred 
points  to  represent  attacks  on  Minuteman  missile  fields  and  a  30  MT 
explosion  at  the  point  marked  by  an  asterisk.  These  events  are  presumed 
to  be  the  total  yield  of  many  single  bombs  detonated  close  together  at 
nearly  the  same  time.  In  Figure  8  we  show  contours  of  the  resulting  dis¬ 
turbances  in  the  same  way  as  was  shown  in  Figure  5.  The  ends  of  the  HF 
link  are  marked  on  the  figures  by  x's.  Distance  on  the  plots  is  distorted 


T  =  15  min 


Figure  5 


Contours  of  the  electron  density  at  250  km  altitude  following 
two  10  MT  explosions.  An  artificial  ionosphere  was  used 
(see  text).  The  tic  marks  are  500  km  apart  and  the  contour 
interval  is  25  km. 
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Figure  5,  (continued) 


Figure  5.  (continued) 
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Figure  6.  Contours  of  electron  density  for  the  event  of  Figure  5 
along  the  line  shown  in  Figure  5.  The  spacing  between 
horizonal  tic  marks  is  500  km.  The  contour  interval  is 


Figure  6.  (continued) 


Figure  7.  An  HF  link  from  Omaha,  Nebraska  to  British  Columbia  and 
the  location  of  explosions  of  100  MT  (star)  and  30  MT 
(asterisk) . 
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T  =  15  min 


Figure  8.  Contours  of  electron  density  (artificial  ionosphere)  at  250  km 
altitude  following  the  attack.  The  X’s  mark  the  end  points  of 
the  link  shown.  The  tic  marks  are  250  km  apart  and  the  contour 
interval  is  50  km. 
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Figure  8.  (continued) 


by  the  necessity  to  project  the  spherical  surface  onto  a  plane.  As  before 
the  artificial  ionosphere  is  used  and  an  altitude  of  250  km  is  shown;  in 
this  case  we  have  shown  the  contours  in  increments  of  50  km  starting  with 
50  km  so  the  E  contour  is  the  ambient  value  at  250  km.  In  Figure  9  we  show 
a  vertical  slice  along  the  line  of  the  HF  link  (as  above,  we  use  the  arti¬ 
ficial  atmosphere  to  clarify  the  results);  contours  are  shown  in  increments 
of  50  km  starting  with  50  km.  In  Figure  10  we  show  the  same  vertical  plane 
as  in  Figure  9  but  in  this  case  we  have  used  a  realistic  ionosphere,  taken 
from  Reference  7,  and  we  have  plotted  contours  of  critical  frequency  start¬ 
ing  with  1  MHz  and  in  increments  of  1  MHz.  This  figure  shows  the  degree  of 

disturbance  in  the  ionosphere  for  the  radio  wave  propagation. 

We  note  in  some  of  these  figures  that  the  model  behaves  in  what 
may  be  an  unphysical  way  when  the  disturbance  becomes  sufficiently  strong. 

To  understand  this  behavior  consider  a  case  where  the  amplitude  of  AZ 
at  its  maximum  value  at  Z  =  300  km  is  greater  than  200  km.  The  model 
would  then  specify  that  the  material  at  300  km  came  from  above  500  km  or 

below  100  km.  But  the  model  specifies  that  the  points  at  100  km  and  500  km 

do  not  move  (that  is,  the  amplitude  of  the  disturbance  goes  to  zero)  thus 
the  model  produces  a  pattern  in  the  ionosphere  which  cannot  be  arrived  at 
by  smooth,  linear,  wavelike  motions.  Close  examination  of  Equation  22 
shows  that  a  problem  of  this  type  will  develop  somewhere  whenever  the 
amplitude  of  AZ  becomes  larger  than  100  km.  The  discontinuity  in  the 
motion  is  near  100  km  and  500  km  where  the  amplitude  of  AZ  is  near  100  km 
and  moves  nearer  300  km  for  larger  values  of  the  maximum  disturbance.  We 
could  remove  the  discontinuity  in  the  underlying  motion  in  various  rela¬ 
tively  simple  ways  but  so  far  we  have  chosen  not  to  do  so.  One  reason 
is  that  the  electron  densities  provided  by  the  model  are  not  discontinuous; 
there  exists  some  motion  which  would  produce  the  predicted  electron 
densities.  Another  reason  is  that  the  problem  occurs  only  for  very  violent 
disturbances  such  that  the  linear  analysis  used  to  develop  the  model  could 
not  be  exp  ;cted  to  give  an  accurate  result.  For  motions  large  enough  to 


Figure  9.  Contours  of  electron  density  along  the  link  (for  the 
artificial  ionosphere)  following  the  attack.  The  tic 
marks  are  250  km  apart  and  the  contour  interval  is 

50  km. 
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Figure  10.  Contours  of  critical  frequency  (realistic  ionosphere) 
along  the  link  following  the  attack.  The  tic  marks 
are  500  km  apart  and  the  contour  interval  is  1  MHz. 
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produce  displacements  in  the  electrons  as  large  as  100  km  the  flow  will 
probably  become  turbulent  and  produce  a  fragmented  ionosphere  at  least  as 
disturbed  as  that  predicted  by  the  model. 

We  have  traced  some  rays  through  the  disturbed  ionosphere  in  this 
model  problem.  These  calculations  were  made  to  demonstrate  that  the  gravity 
wave  disturbance  produces  considerable  effects  on  the  propagation  of  HF 
radio  waves.  We  have  not  yet  made  calculations  in  enough  detail  to  analyze 
the  expected  performance  of  a  system  on  the  proposed  link.  By  tracing  rays 
in  the  undisturbed  ionosphere  we  found  that  a  25  MHz  ray  launched  at  about 
6°  elevation  would  propagate  to  about  2800  km  on  a  single  hop.  15  MHz  rays 
would  not  propagate  to  2800  km  on  a  single  hop  but  for  an  elevation  angle 
of  about  18°  two  hops  propagated  to  2800  km.  We  then  ran  a  few  rays  in  the 
very  disturbed  environment  at  60  minutes.  The  15  MHz  rays  launched  near  18° 
elevation  all  came  to  earth  behind  and  to  the  left  of  the  transmitter.  We 
made  no  further  calculations  at  15  MHz  at  60  minutes.  The  25  MHz  ray  launched 
at  6°  elevation  went  only  1700  km,  and  no  ray  launched  along  the  line  of  the 
link  went  more  than  2765  km.  Rays  launched  near  1°  elevation  and  somewhat  to 
the  right  of  the  intended  receiver  did  land  in  the  region  around  2800  km  along 
the  link.  Whether  this  is  the  only  launch  angle  to  arrive  at  this  point  we  do 
not  know,  but  we  did  notice  an  interesting  phenomena  at  about  2400  km.  In 
Figure  11  we  show  the  landing  point,  of  two  different  bundles  of  rays.  Fig¬ 
ure  11a  shows  a  bundle  centered  at  a  launch  elevation  of  2.70°  and  a  launch 
azimuth  of  325.50.  In  Figure  lib  we  show  the  bundle  centered  at  a  launch 
elevation  angle  of  7.93°  and  a  launch  azimuth  angle  of  324.90°.  In  each  plot 
the  x  in  the  center  is  the  same  point  some  2400  km  along  the  path  shown  by 
the  line  in  Figure  7.  As  some  sort  of  comparison  we  show,  in  Figure  11c, 
the  landing  points  of  the  bundle  of  rays  at  the  lower  launch  elevation  angle 
for  the  ambient  ionosphere.  The  center  of  the  plot  must,  of  course,  be  dif¬ 
ferent  but  a  comparison  of  the  patterns  of  Figure  11c  with  those  of  Figures  11a 
and  lib  shows  the  extreme  distortion  in  the  ray  paths  produced  by  the  gravity 
wave.  While  we  have  not  been  able  to  analyze  the  signal  at  the  location  of  the 
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center  of  Figures  11a  and  lib  in  this  report, 
path  problem  may  be  present  for  some  systems. 


it  appears  likely  that  a  multi 


We  now  consider  the  time  of  240  minutes  as  shown  in  Figure  10  and 
again  consider  the  case  of  15  MHz  and  two  hops.  We  have  calculated  ray  path 
for  rays  launched  along  the  line  in  Figure  1  and  having  elevation  angles 
from  zero  to  twenty-five  degrees.  The  rays  all  land  fairly  near  the  dir¬ 
ection  of  launch,  up  to  a  few  degrees  deviation,  so  we  can  give  some  feel 
for  the  results  by  plotting  the  range  vs.  the  launch  angle;  this  has  been 
done  in  Figure  12.  The  transmitted  energy  is  seen  to  become  highly  con¬ 
centrated  in  the  vicinity  of  3400  km.  It  appears  that  even  after  4  hours 
the  gravity  wave  is  causing  interesting  radio  wave  phenomena  at  HF  frequen¬ 
cies.  It  should  be  possible  to  perform  an  analysis  of  these  effects  from 
a  systems  point  of  view  in  the  near  future. 
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Figure  12 


The  range  to  which  rays  launched  at  elevation  angle  6 
propagate  for  15  MHz,  two  hops,  for  the  240  minute 
envi ronment. 


59 


SECTION  5 
CONCLUSIONS 


In  this  report  we  have  described  a  model  which  gives  the  ionos¬ 
pheric  effects  of  the  gravity  waves  produced  by  low  altitude  nuclear  ex¬ 
plosions.  For  large  distances  from  the  burst  the  model  is  based  largely 
on  data;  for  the  region  near  the  burst  the  model  is  based  entirely  on 
calculations.  The  model  is  described  by  a  single  set  of  equations,  how¬ 
ever,  and  is  thus  continuous  with  no  "switch  points"  which  might  lead  to 
anomalous  behavior.  The  gravity  wave  model  has  been  programmed  into  a 
ray  trace  computer  code  giving  us  the  ability  to  analyze  the  changes  the 
disturbed  ionosphere  produces  on  HF  radio  waves. 

The  gravity  wave  effects  can  be  quite  severe  in  the  case  of  a 
large  nuclear  attack  and  strong  effects  last  for  periods  of  several  hours. 

In  the  present  work  we  have  been  able  to  do  only  a  few  example  calculations. 
These  examples  confirm  that  the  gravity  wave  model  produces  significant 
ionospheric  distortions,  as  expected,  that  the  ray  trace  codes  can  calculate 
the  raypaths  through  the  distorted  ionosphere,  and  that  the  gravity  wave 
phenomenon  severely  distorts  the  raypaths  for  times  up  to  at  least  several 
hours  after  a  nuclear  attack. 
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APPENDIX  A 

(Taken  from  DNA  3494T) 
SHOCK  WAVE  IN  THE  ATMOSPHERE 


DECAY  OF  THE  BLAST  WAVE 


After  a  nuclear  explosion  a  substantial  portion  of  the  yield, 
perhaps  30%,  is  lost  from  the  region  immediately  around  the  burst  by  radiation 
and,  except  for  some  heating  in  the  air,  plays  no  further  role  in  the  hydro- 
dynamic  effects  which  follow.  The  remaining  energy  is  chiefly  internal  energy 
in  the  high  temperature  fireball  which  is  approximately  an  isothermal  sphere. 
As  this  sphere  rapidly  expands  it  drives  a  very  strong  blast  wave  ahead  of 
it.  As  a  consequence  the  original  internal  energy  of  the  small  hot  region 
is  converted  partly  to  the  energy  of  the  blast  wave,  while  the  remainder  is 
in  the  internal  energy  of  the  fireball  which  has  expanded  to  pressure  equi¬ 
librium.  To  get  an  estimate  of  the  fractions  involved  consider  a  small  hot 
region  of  ideal  gas  of  density  p  and  pressure  P  which  adiabatically 
expands  to  a  new  state  p",  P'  with  adiabatic  exponent  y.  If  the  region 
contains  mass  M  ,  then  the  original  and  final  internal  energies  are 


E  =  M 


P 

(y-l)p 


and  E  "  =  M  •  — 7-r — *- 
(Y-1)P 


(A-l) 


By  the  adiabatic  condition 


(P/P)  P 


(1-Y) 


(P7p')(p')(1'Y) 


C  A- 2  ) 


so 


E7E  =  (p7p)y_1 


(A-3) 


A-l 


A  typical  case  of  interest  might  be  y=1.3  and  pVp  ~  1/20  which  gives 
E'/E  »  1/3.  That  is,  2/3  of  the  original  internal  energy  would  be  transferred 
to  the  blast  wave.  Consequently,  perhaps  50%  of  the  yield  resides  in  the 
strong  blast  wave.  This  energy  is  deposited  in  the  region  through  which  it 
passes  and  it  is  of  considerable  interest  to  know  just  where  and  how  this 
deposition  occurs.  In  order  to  shed  some  light  on  this  question,  let  us 
consider  first  the  more  accessible  question  of  how  the  pressure  of  the  blast 
wave  declines  with  distance. 


Pressure  Decline 


For  the  case  of  a  point  explosion  there  exists  an  analytical 
solution  for  the  complete  flow  field  if  the  external  pressure  is  negligible 
compared  with  the  shock  pressure.  The  pressure  at  the  shock1  is  given  by 


D  -  P _  JlQ.  _  D~ 

s  '  25 (Y+l)  K(y) 


CA-4) 


where 


with 


13y2-7y-*-12 

5(2y+1)(3y-1) 


(A-5) 


A-2 


which  is  the  order  of  1  for  y  =  1.3.  E0  is  the  energy.  Since  P«R  3,  it 
is  apparent  that  the  sphere  functions  as  a  unit  and  the  shock  will  not  be 
able  to  run  away  and  leave  a  detached  inner  region.  In  a  realistic  case 
the  strong  shock  condition  must  eventually  cease  to  be  satisfied  due  to  the 
finite  pressure  Po  of  ...  >e  ambient  medium.  Still  later  when  the  shock 
weakens  to  the  acoustic  level  the  result  should  become  the  usual  acoustic 
form 


Ps  =  P0  +  a/R.  (A- 6) 

The  transition  between  these  two  forms  is  of  course  described  by  hydrodynamics 
although  no  analytic  solution  exists.  Brode2  has  made  a  study  of  this  question 
by  use  of  a  finite  difference  numerical  method  and  the  resulting  pressure 
versus  range  curve  is  shown  in  Figure  1.  The  initial  values  for  Brode' s 
calculation  were  taken  from  the  analytic  point  source  solution  with  a 
dimensionless  length  unit  of 


where  Rg  is  the  shock  radius,  E-j-qt  t*ie  and  p°  t^e  ambient  pressure. 

As  expected  there  is  a  R  3  behavior  until  the  shock  pressure  ceases  to  be 
much  larger  than  ambient. 

Measurements  of  overpressure  versus  range  have  been  made  at  most 
of  the  U.S.  nuclear  tests  .  In  order  to  compare  such  results  with  calculations 
it  is  necessary  to  suitably  scale  the  data  taking  into  account  the  yield, 
altitude  and  physical  surroundings  of  the  bursts.  This  will  be  done  in  the 
usual  way.  In  particular,  suppose  a  measurement  was  made  from  a  burst  with 
yield  Y,  ambient  pressure  P0  ,  density  p0,  and  sound  speed  c0,  and 
determined  a  pressure  P,  density  p,  or  velocity  u,  radius  R  at  time 
t.  Then  the  scaling  laws  indicate  that  values 


wwsm i  aw-  -■  ••  *  - 


A-3 


will  be  found  at 


and 


where 


(A-8) 


(A-9) 


for  conditions  indicated  by  a  new  set  of  parameters  indicated  by  primes. 

Here  e  is  a  parameter  indicating  the  surroundings  of  the  burst 

2  burst  over  water 

1.5  burst  over  land  (A- 10) 

1  free  air  burst 

Data  from  a  large  number  of  US  tests,  scaled  to  1  kiloton  TNT  equivalent 
(1  KT  E  4.185xl019  ergs),  are  shown  on  Figure  2.  The  results  of  the  finite 
difference  calculation  appropriately  scaled  (ET0TAL  =  1  KT)  are  also  shown 
and  the  agreement  is  excellent.  Much  of  the  data  for  low  overpressures  is 
not  very  good  in  that  it  comes  from  high  yield  bursts  and  .-o  corresponds  to 
measurements  taken  near  the  ground  after  the  blast  wave  has  traveled  a  very 
great  distance.  This  would  introduce  ground  effects  and  variations  with  local 
conditions.  For  example  AP/P0  =  .05  for  a  yield  of  1  Megaton  (MT)  at  sea 
level  occurs  at  a  distance  of  about  20  kilometers. 
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•  Data  scaled  to  1  KT 

°  Low  yield  data  scaled  to  1  KT 

*  HA  data  scaled  to  1  KT  i 

®  Brode's  results  scaled  to  1  KT 
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Figure  A-2.  Pressure  versus  Range  for  Test  Data  Scaled  to  1  KT. 


Acoustic  theory  predicts  that  AP/P0  ~  R~  1  for  large  R.  However 
this  corresponds  to  assuming  that  no  energy  is  lost  from  the  pulse  which  is 
not  true.  Bethe  has  developed  a  semi-acoustic  theory  in  which  lowest  order 
non-linear  corrections  are  introduced  in  the  acoustic  equations  and  the  lowest 
order  entropy  gain  term  is  taken  from  the  Hugoniot  shock  conditions.  Assuming 
that  the  shape  of  the  pulse  is  constant,  Bethe1  was  able  to  show  that  for 
spherical  symmetry  the  asymptotic  form  is 

AP/P0  ~  (R  An(R/R0))_1  .  (A-ll) 

Such  behavior  is  shown  in  Figure  1  and  it  appears  that  this  asymptotic  behavior 
dominates  for  overpressures  less  than  roughly  .05. 

As  the  blast  wave  moves  out  away  from  the  region  approximated  by 
a  strong  shock  from  a  point  source,  the  shape  of  the  pulse  becomes  that  of 
Figure  3.  The  sequence  of  changes  leading  to  that  shape  are  shown  in 
Figure  4.  Clearly  the  limiting  form  can  occur  only  when  pulse  has  an  over¬ 
pressure  less  than  1.  In  Figure  5  the  speeds  corresponding  to  the  pressures 
of  Figure  4  are  shown.  The  particle  speed  in  the  pulse  behaves  much  as  the 
pressure.  It  is  to  be  expected  that  the  shapes  of  u  and  P  versus  r 
would  be  similar  in  that  for  acoustic  theory  we  know 

u  =  (P-P„)/p0c0  (A- 12) 

where  Po  is  the  ambient  density  and  c0  the  sound  speed. 

There  exists  a  large  quantity  of  data  on  the  positive  pulse 
duration  but  much  of  it  is  of  poor  quality.  The  data  scaled  to  1  KT  at  sea 
level  are  shown  on  Figure  6.  The  scatter  of  the  data  at  the  larger  radii 
shown  is  clearly  indicative  of  the  existance  of  extraneous  effects  for 
propagation  over  long  paths.  Special  attention  should  be  paid  to  the  points 
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Relative  Radius 

Figure  A-4.  Development  of  Pressure  Pulse  in  Brode’s  Solution.  The 
Variables  are  Scaled  to  their  Values  at  the  Shock. 


Radius  (meters) 

Figure  A-6.  Positive  Pulse  Duration  versus  Radius  Scaled  to  1  KT  for 
Data  and  Fit  as  given  in  Text. 


from  shot  HA  for  which  the  measurements  were  taken  in  free  air  by  attaching 
instrument  packages  to  parachutes  which  were  released  near  the  burst  altitude. 
Note  that  the  points  on  Figure  2  for  the  pressure  of  HA  all  lie  below  the  fit 
curve.  If  the  assumed  yield  was  arbitrarily  increased  from  3  KT  by  about 
70%  the  pressure  points  would  be  brought  onto  the  fit  curve.  The  same  change 
would  bring  the  points  for  HA  on  Figure  6  up  about  20%  and  provide  a 
reasonably  consistent  set  of  data.  The  line  drawn  to  provide  a  fit  to  the 
data  has  been  chosen  to  be  of  the  form 


D+(R)  =  Do 


co  2y 


1+  7T~  tn(R/R0)* 


(A- 13) 


Dn 


which  is  based  on  conditions  which  are  discussed  in  the  section,  "Model  for 
Propagation" , 


Energy  Decline 


The  energy  in  the  outward  traveling  pulse  is 


E  =  4n  J  R*dR  ~  f  RZdRPuZ  ,  (A- 14) 

which  is  the  sum  of  the  internal  energy  (IE)  and  kinetic  energy  (KE) .  Since 
for  the  acoustic  case 

P-P0  =  upc  (A- 15) 

it  might  be  supposed  that  the  KE  is  of  lower  order  than  the  IE.  However, 
this  is  not  so  due  to  the  negative  phase  of  the  pulse.  In  fact  for  the 
acoustic  case  KE  =  IE.  For  the  case  of  a  strong  shock  it  continues  to  be 
the  case  that  KE  =  IE  at  the  front.  The  Rankine-Hugoniot  conditions  expressing 
mass,  momentum,  and  energy  conservation  across  the  shock  require  that  for 
a  strong  shock. 


A-12 


(A- 16) 


Ps/(Y-U  -i  PsuJ  . 

Consequently  we  will  consider  that  E  =  2xKE  for  all  our  cases. 

The  energy  in  the  pulse  shown  in  Figure  3  may  be  approximated  as 
that  in  the  similar  linear  form.  Then 


ES  p4ttR*  (L_)(u')a  +  4ttR2  (L+)  (u+)2p 

=  E'  +  E+  (A-17) 

Relative  dimensions  for  the  pulses  in  the  well  developed  negative  phase 
region  show  that 


and 


P+  2  3P_  or  u+  —  3u 


L  S  2.5  L 

+ 


(A- 18) 


so  E+/E  -  .78.  Consequently  we  will  hereafter  consider  only  the  energy  in 
the  positive  phase  as  the  energy  of  the  pulse  and  the  general  expression 
used  will  be 


5+  =  477  J 


R2puzdR  . 


R  -L 
s  + 


(A-19) 


For  the  purpose  of  computing  E+  from  data  we  will  presume  that  the 
pulse  is  linear  and  weak.  Then 


A- 1 3 


and 


AP  =  up0c0  , 


PSPo  so 


CA-20) 


The  length  of  the  pulse  L+  is  related  to  the  passage  time  of  the  pulse  At 
by  L+=c0At  since  the  shock  is  weak  and  traveling  at  nearly  the  sound  of 
speed  Co  .  Figure  7  shows  the  widely  scattered  data  scaled  to  1  KT  at  sea 
level.  Three  sets  of  data  taken  from  the  same  tests  are  shown,  one  being 
from  shot  HA  with  the  yield  increased  by  70%  as  discussed  previously.  These 
sets  provide  a  better  feeling  for  the  radial  behavior  of  E+ZEjoTAL  when 
taken  individually. 


It  is  apparent  that  by  the  time  the  shock  reaches  an  overpressure 
of  about  .05  only  5%  of  the  yield  resides  in  the  positive  pulse  part  of  the 
shock.  In  order  to  understand  this,  the  mechanism  by  which  energy  is  lost 
must  be  investigated. 


The  sequence  of  events  which  leads  to  the  loss  of  energy  from  the 
pulse  is  as  follows: 

1.  Material  is  shocked  increasing  entropy,  density,  pressure, 
and  temperature. 

2.  After  passage  of  the  shock  the  material  reexpands,  approximately 
adiabatically. 

3.  The  material  finally  reaches  the  original  pressure  and  ceases 
to  expand  but  it  is  now  at  higher  than  original  temperature. 

4.  Since  the  pressure  is  ambient  no  further  expansion  occurs  and 
the  increased  internal  energy  is  no  longer  available  to  the 
shock . 


This  will  be  referred  to  as  the  Taylor3  loss  mechanism. 
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Scaled  data 

Scaled  data  from  one  burst 
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Figure  A-7.  Percent  Yield  in  Positive  Pulse  versus  Radius  Scaled  to 
for  Data  as  Compared  to  Taylor  Loss  Result. 


The  increase  in  internal  energy  per  unit  mass  will  be 


(Y-l)Pf  (Y-l)Po  (Y-l)Po 


(po/p£  -  l) 


(A- 21} 


where  is  the  final  density  after  reexpansion  which  is 


pf  =  ps 


(A-22) 


The  subscript  s  denotes  values  at  the  shock  position.  Thus  by  the  Rankine- 
Hugoniot  relations 


Ps/Po  = 


Y" 1  +  2/M2 


Vp°  -  1  +  (m2-^ 


(A-23) 


where  M  =  V^/c,,  is  the  Mach  number  of  the  shock.  The  entropy  increase 
per  unit  mass  due  to  the  shock  will  be 


AS=  ^[*n(psps-Y)_£n(poPf-Y) 


(A-24) 


where  P  =  PR*T.  Therefore  the  energy  per  unit  mass  which  can  no  longer 
be  available  for  work  will  be 


ter  m) 


(A- 25) 


(Y-l)Po  ln  ^o/Pf) 
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Since  Po/p^  ■*  1  for  weak  shock  and  £n(l+e)-*-e  one  sees  that 
AEs  =  yAI  for  AP/P0«1.  If  one  takes  this  limit  it  is  found  that 


Al 


(Y+1) 

12  y3 


(A -2b) 


Note  that  the  first  two  orders  do  not  appear.  On  Figure  8  yAIo>  yAI  , 

AI,  and  Es  are  plotted.  The  energy  which  is  lost  from  the  shock  will  be 

at  least  equal  to  the  larger  of  AEg  and  AI .  For  shocks  that  are  not  too 
AP 

strong  say  ^—  <  10,  a  suitable  expression  for  the  rate  at  which  energy  is 
“  o 

lost  from  the  shock  is 


_1 

4’r2^  p.((^)(fc)Y  -  ■)  •  <A-27> 

From  this  one  can  find  the  shape  of  Egfloc|c  versus  R  curve  given 
Ps(RJ.  However,  there  is  no  firm  initial  condition  with  which  to  normalize 
this  curve.  On  Figure  7  a  set  of  Es^qc^(R)  curves  are  given  based  on 
Ps(R)  from  Figure  2.  It  appears  that  there  is  reasonable  agreement  with 
the  data  although  the  data  are  poor. 


dE  ,  , 

shock 

dR 


If  one  assumes  a  point  source  solution  it  is  possible  to  find 
accurately  the  energy  which  is  available  to  the  shock.  Using  Equation  (2-27) 
in  conjunction  with  Ps(R)  from  the  point  source  solution  from  Equation  (2-4), 
one  finds  that,  E^,  =  .23  E^  . 


P  =  108. 
s 

Y  =  1.2 
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Note  that  y=  1.2  has  been  used  due  to  the  high  temperature  influence  on  the 
equation  of  state.  Such  a  calculation  cannot  be  taken  too  seriously  because 
the  point  source  solution  is  not  applicable  to  a  real  explosion  at  early 
times  due  to  the  radiation  phase  of  real  explosions. 

Still,  a  value  of  .23,  which  would  be  at  about  85  meters  on  Figure  7, 
is  in  the  range  of  the  data.  Perhaps  the  energy  loss  by  radiation  in  the 
real  case  is  of  the  same  magnitude  as  the  Taylor  loss  energy  for  the  very 
strong  shock  at  early  times  in  the  point  source  model.  Note  that  by  the 
point  source  model  ^sjloci.  =  >64  e-j-qtal  whereas  for  the  realistic 

P  =  1000  P0 
s 

Y  =1.2 

case  the  radiation  phase  would  still  be  important  allowing  energy  to  be 
easily  transported  throughout.  Thus  some  energy  which  would  have  been  lost 
from  the  shock  can  return  to  the  shock  by  radiation. 

s 

In  conclusion  it  seems  certain  that  the  portion  of  the  yield  which 
is  in  the  horizontally  outgoing  weak  shock  wave  cannot  be  more  than  5%  when 
the  pressure  has  fallen  to  5%  of  ambient.  For  the  region  beyond  AP/P0  =  .05 
the  pressure  Ps(R)  appears  well  described  by  the  semiacoustic  result 

P_(R)-P0  =  ■  ~  •  (A-28) 

R/£n(R/R0) 

As  a  consequence  the  energy  content  decreases  like  [£n(R/R0)]  1  at  large 
distances. 


BEHAVIOR  IN  A  STRATIFIED  MEDIUM 


Previous  discussions  were  concerned  with  propagation  through  an 
isotropic  mediivr.  However,  for  propagation  in  the  atmosphere  this  cannot 
be  a  suitable  description  unless  the  characteristic  length  (E^,q^,/Po)  /3  <<  H 
where  H  is  the  scale  height.  For  cases  of  current  interest  this  condition 
is  not  satisfied.  Consequently  we  must  consider  the  effect  of  the  stratified 
atmosphere  on  propagation  of  a  shock  wave. 

Acoustic  Propagation 


In  order  to  begin  to  understand  such  effects,  let  us  consider  the 
propagation  of  an  initially  spherical  pulse  into  an  isothermal  atmosphere 


which  is  stable  under  a  constant  gravitational  force  acting  in  one  Cartesian 
dimension  as  described  in  the  linearized  (acoustic)  equations  of  hydrodynamics. 

In  general. 


W  +  ^  =  •  I  VP'gSz  ' 
~  +  V*Vp  =  -  pV*V  ,  and 

3P  -+ 

~  +  V* VP  =  -  yPV* V  , 


(A-29) 


where  the  third  equation  is  equivalent  to  the  adiabatic  condition 


_d_ 

dt 


(pp't)-  0 


(A-30) 


=>  ■.  'v  -  -VJrl**  rT  Ifti 
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Now  for  an  ideal  gas 


P0  =PoR*T 

so  that  the  steady  state  unperturbed  condition  is 

p  -  p  e_  ) 

sL  6  / 

-  z/H  I 
Po  =  PsL  6  ) 

where  H  =  -  —  =  ~  —  =  -  R*T0  . 
g  Po  g  PsL  g 

The  linearized  equations  then  give 
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cos0  „  ,r  sin0  _ 
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Po 
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(A-32) 


(A-33) 
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where 


V  =  VT&T  *  V0ee  ,  p=  p0  +  P,  and  Pi  =  P0  +  Pi  • 


If  we  temporarily  assume  Vg  =  0  and  write  =  U/r 


for  U  is  found  as 


32U  2  32U  2  cos6  3U 


at 


c 

2  0 


3r 


2  C  0  H  3r 


Therefore  for  a  spherical  outgoing  wave  with 
U  -  g1 (kr-wt) 

we  find  the  dispersion  relation 


k2  +  i 

c'i  n 


or 


or 


,  -i  cos0 
k  =  —  -  —  + 


c 


2H 


1  - 


A  f  4a)2  cqs29  I 

2  LT'  H2  J 

] 


c2  cosz0 1 

_ 0 _  I  ;  cos 

4HZ  w2 


The  limit  H  -*■  °>  gives  the  usual  isotropic  case.  However  if 


W  <  2H  COS0 


k  is  purely  imaginary  and  such  frequencies  will  not  propagate, 
reason  (jJc  =  ^  is  known  as  the  acoustic  cutoff  frequency.  The 
value  in  the  lower  atmosphere  is  -2*10  2  sec  1  . 


an  equation 

(A-34) 


(A-3S) 


(A-36) 


(A-37) 

For  this 
numerical 
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The  velocity  of  the  particles  will  therefore  be 


V 

r 


<*  —  exp  i  ((Re  k  )r-tot)  exp 


I  COS0 
\  2H 


') 


(A-38) 


Thus  for  an  upwardly  traveling  wave  the  particle  velocity,  and  also  the 
pressure,  increases  exponentially  according  to 

r 


V  ^  Co  —  =  — -  e 
r  0  p0  c0po  r 


2H 


(A-39) 


in  the  acoustic  limit  for  to  » to 

c 


Returning  to  the  question  of  VQ  ,  we  find  using  our  solutions. 


(A-40) 


Therefore  if  —  »^r  it  will  be  true  that  |VQ|«|V  I  and  the  V„  =  0 
c  2H  1  0 '  1  r 1  0 

assumption  (strictly  radial  flow)  is  justified.  The  condition  for  validity 
is  perhaps  better  written 


A  «  47iH  (A -41) 

where  A  =  2?rc/co  is  the  wavelength.  Note  that  the  condition  for  upward 
propagation  (2-37)  could  have  been  written  — <  4ttH  although  since 

A  =  2ir/Re  k  =  2ttc0/w  [l-co^/co2]'*5  ,  (A-42) 

the  quantity  2ttc0/w  ceases  to  be  the  wavelength  unless  oo»u)c. 
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The  important  conclusion  is  that  so  long  as  oj»ojc  ,  which  is 

usually  the  case,  the  pressure  of  a  pulse  propagating  up  into  a  decreasing 

1  r /2H 

stratified  medium  will  vary  as  Pi/Po'*  —  e  x  according  to  the  acoustic 
equations.  Of  course  no  energy  will  be  lost  from  the  pulse  and  the  pulse 
will  not  change  shape  so  long  as  oo  >>  coc . 

Non-Linear  Propagation 

There  exists  no  analytic  solution  to  the  general  non-linear  problem 
of  propagation  of  a  pulse  into  a  medium  of  variable  density.  However,  some 
related  problems  have  been  studied  which  contain  some  aspects  of  the  problem 
of  current  interest4.  Raizer5  produced  a  similarity  solution  for  the  case 
of  a  strong  plane  shock  in  an  infinite  medium  with  exponentially  varying 
density.  This  solution,  also  found  by  Grover  and  Hardy6,  allows  the  shock 
wave  to  reach  infinity  in  a  finite  time  with  a  shock  velocity  at  position  x 

Vs(x)  =  c  eaX/H  (A-43) 

where  a  is  a  function  of  y  given  in  Figure  9.  Grover  and  Hardy  also 
studied  the  problem  by  a  finite  difference  code  and  found  that  asymptotically 
the  solution  approached  the  similarity  solution  largely  independently  of  the 
initial  pulse  configuration.  The  similarity  character  of  this  solution  is 
such  that  the  ratio  of  the  variables  p,  P,  and  u  behind  the  shock  to  that 
at  the  shock  is  constant  for  a  fixed  distance  behind  the  shock.  Consequently 
one  cannot  speak  about  the  length  of  the  pulse. 

They  considered  by  their  finite  difference  code  the  same  problem 
in  spherical  and  cylindrical  geometry  and  found  the  asymptotic  behavior  to 
involve  the  same  exponential  factor.  This  is  shown  by  Figure  10. 
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ntial  Growth  Factor  for 
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■.  '***- 


Shock  Velocity  (units  of  C) 


Figure  A-10.  Shock  Velocity  versus  Position  for  y  =  2  (from  Grover  and  Hardy). 
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In  a  similar  vein  Enstrom  and  Brode8  have  considered  finite 
difference  solutions  for  both  weak  and  strong  shocks  in  planar  geometry. 
For  strong  shocks  the  growth  is  very  much  like  that  of  Grover  and  Hardy7, 
that  is 


V 

s 


oc  e 


ax/H 


(A-44) 


while  for  weak  shocks  the  growth  is  much  more  rapid  as  indicated  on  Figure  9. 
Strong  and  weak  here  mean  whether  P*H  is  much  less  than  or  much  greater 
than  the  energy  per  unit  area  of  the  shock. 


From  these  considerations  we  expect  that  the  upwardly  traveling 
shock  from  a  nuclear  explosion  will  behave  like 


V  «  — 
s  r 


3a 

+  -  T 

.  H  r 


(A-45) 


once  the  shock  wave  gets  sufficiently  far  away  from  the  burst  point.  Since 
H  is  a  function  of  height  and  increases  rapidly  above  100  km,  one  may 
expect  the  growth  to  decrease  substantially  in  the  real  atmosphere  in  this 
region. 


MODEL  FOR  PROPAGATION 

In  order  to  understand  how  the  upwardly  propagating  shock  varies 
with  yield  and  burst  altitude,  it  would  be  useful  to  have  a  simple  model  of 
the  process  short  of  the  use  of  expensive  two  dimensional  finite  difference 
codes.  Otterman8  has  suggested  a  model  for  such  propagation  of  nearly 
acoustic  pulses  and  Reed 9  has  provided  an  additional  condition  which  allows 
some  loss  of  energy  as  in  the  semi -acoustic  theory  of  Bethe1 . 

The  Otterman-Reed  model  provides  for  a  triangular  pressure  pulse 
propagating  with  length  L(R)  and  relative  overpressure  3(R)-  This  pulse 


A-27 


is  to  be  such  that  the  flow  is  strictly  radial  so  that  the  model  is  one 
dimensional  and  spherical.  It  is  assumed  that  the  tail  of  the  pulse  which 
is  at  ambient  pressure,  travels  at  the  ambient  sound  speed  while  the  head 
of  the  pulse  advances  at  the  shock  speed  ,  and  that  the  shape  of  the  pulse 
remains  triangular.  This  means  that 


^=Vs(R)-c0  •  (A-46) 

It  is  reasonable  that  the  tail  advances  at  the  sound  speed  since  the  tail 
is  at  the  ambient  condition  where  c  is  exactly  the  speed  with  which 
disturbances  are  propagated.  However,  the  assumption  about  constancy  of 
shape  requires  some  support.  The  energy  in  the  pulse  is  taken  as 

Es  =  ^R2B2LPo(R)  (A-47) 

which  is  obtained  for  an  acoustic  spherical  pulse  with  a  triangular  pressure 
distribution  of  length  L.  This  is  readily  found  by  taking  Eg  as  twice  the 
kinetic  energy 


Es  =  4ttR  Pq  |  uzdR 


"Jui 


(A-48) 


where  B  =  uPoc  .  To  complete  ‘he  model  Reed  observed  that 


dis  =  .  IHD  p  g2R2  =  .  Hill  I  E 

dR  3  Yz  0  4y  L  s 


(A-49) 


for  the  small  0  limit  as  noted  earlier.  Finally  one  notes  that  the 
Rankine-Hugoniot  condition  on  the  shock  velocity  in  the  weak  shock  limit  is 
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Therefore  the  final  set  of  equations  is 


and 


dL  _  Y+l  R 
dR  ~  4y 

1  dEs  Y+l  3 

E  dR  =  ‘  4y  L 
s 


(A-51) 


where  E  =  R232LP0 (R) .  Taking  D+(R),  the  duration  of  the  positive  pulse, 

S  3Y  +  + 

to  be  given  by  D  (R)  =L  /c0,  the  solution  to  the  above  equations  gives  the 
result  previously  quoted  in  Equation  (2-13) . 


This  model  is  not  suitable  for  our  purposes  because  it  is  limited 
to  very  weak  shocks.  However,  the  basic  ingredients  are  not  explicitly 
dependent  upon  the  weak  shock  assumption  so  one  can  attempt  to  modify  the 
model  for  higher  pressure.  Kahalas  and  Murphy10  have  added  the  second  order 
corrections  to  the  shock  velocity  and  energy  loss  which  allow  a  somewhat 
stronger  shock  to  be  described,  but  this  is  still  inadequate  for  our  purposes. 
We  propose  to  avoid  expanding  in  powers  of  B  and  simply  express  the  shock 
speed  and  energy  loss  in  their  general  form.  That  is  we  will  take  as  our 
model  equations 


and 


here  pg  and  Vg  are  found  from  the  Hugoniot  conditions  as 


Ps/Po  - 


Y+l 


Y-l  + 
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(A-53) 


The  quantity  dE^/dR  will  represent  energy  lost  to  the  gravitational 
potential . 


The  energy  in  the  pulse  will  be  taken  to  be  twice  the  kinetic 
energy.  It  will  be  assumed  that  the  pulse  is  triangular  in  the  particle 
velocity  (rather  than  pressure)  as  viewed  in  Lagrangian  coordinates.  Thus 
the  energy  will  be 


Po (r)  u  r2dr 


(A-S4) 


where  u  = ug  j  and  p0(r)  = p*exp(-r/H) .  P*  is  the  ambient  density 

at  the  burst  altitude.  Generally  H(r)  is  a  function  of  position  but  it 
is  only  slowly  varying.  Note  that  H  is  the  average  scale  height  up  to  r, 
not  the  local  scale  height.  Assuming  H  can  be  removed  from  the  integral 
we  define  a  function  f  by 


T1  (ff 


LRZP0(R)  exp(fL/H) . 


(A-55) 


Thus  f  is  the  fraction  of  L  back  from  the  pulse  front  at  R  such  that 


R  K 

J  r*  e-r/H  (r^)2  dr  .  „z  e-(R-fL)/Hj 
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- - -  ---y — Ki>M  - 


This  requires 


where 


x  = 
h  = 
g  = 


L/H,  y  =  R/H, 

24  + 12(y-x)  +  2(y-x)2  ,  and 

e  x  (  x1*  +  4x3  +  12x2  +  24x  +  24  +  2(y-x)  (x3  +  3x2  +  6x  +  6) 
I  +  (y-x)2(x2  +  2x  +  2) 


(A-57) 


J 


The  particle  velocity  at  the  shock  us  is  related  to  g  by  the 
Hugoniot  conditions  giving 


S  p 

T  =  y/1  +  g(y+l)/2y  1  (A-58) 

Note  that  for  (3->-0  ,  ug  «  $  and  for  g-*00,  ug  “/g- .  It  may  appear  that 
us  <*  g  is  inconsistent  with  the  assumption  of  the  equality  of  the  internal 
and  kinetic  energies.  In  fact  this  demonstrates  the  existence  of  the 
negative  phase  for  weak  pulses  which  gives  a  net  contribution  to  the  intfemal 
energy  which  is  proportional  to  g2. 


The  loss  of  energy  to  gravity  will  be  found  by  noting  that  following 
the  passage  of  the  shock  the  atmosphere  will  return  to  ambient  pressure  P0 
but  to  a  temperature  T>T0  .  For  an  ideal  gas  this  will  be 


(A-59) 
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with  a  density 


Pf  *  P0  -f 


(A-60) 


Taking  R(r,t)  as  the  Eulerian  coordinate  the  continuity  condition 


is 


pfR2dR  =pQ (r)  r2dr 
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so  that 


R3  (r)  -r3  = 
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r2dr 
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The  amount  by  which  the  atmosphere  originally  at  r  will  be  raised  is  thus 


Ar  =  R(r)  -  r 


and  the  energy  lost  from  the  pulse  will  be 


dE 


jf-  =  4irrz  p0(r)  (R(r)-r)  g 
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From  Equation  (2-55)  following  differentiation  with  respect  to  r  we  obtain 
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N°w  ~  is  known  from  the  atmosphere  and  — -  is  known  from 

Equation  (2-58) .  The  function  f  is  readily  differentiated  to  give 
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Combining  these  results  we  find 


4ity 
'  Y-l 

4%)  te) 

1 

Y  . 

\  dE  ' 

f 

) 

/!  +  f  +  _ L 

3f 

\+  dB  / 

2 

= 

dr 

U  H  H2 

3x 

1  dr  \ 

(us/c) 

+  1  - 

1  +  r  dH  L 

/ 

L_  dH  _f 

+  M  / 

r 

H  H  dr  H 

l 

dr  x 

9y  \ 

fL 

dH 

‘  H2 

dr 

d(us/c) 


(A-66) 


All  parameters  are  either  known  or  functions  of  8  and  L.  Therefore  if 
we  use  this  in  conjunction  with 

dL  ‘VC)-1 


dr  =  (V  /c)  (A-67) 

we  have  a  pair  of  integrodifferential  equations  for  B(r)  and  L(r).  Note 
that  dEg/dr  depends  on  the  history  of  0. 

SOLUTION  OF  THE  MODEL  EQUATIONS 

These  equations  are  readily  integrated  numerically  using  a  Runge- 
Kutta  procedure  for  initial  values  0o(ro)  and  L  (ro)  provided  dE^/dr 
is  suitably  handled.  The  results  are  relatively  insensitive  to  dE^/dr  so 


it  is  treated  by  taking  Ar  to  be  constant  in  each  interval  of  a  scale  height 
which  is  taken  from  the  integration  through  the  previous  scale  height.  That 
is  we  use 
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for  ro  +  (n-l)Ax  < r < r0 + nAx  where  Ax  has  been  chosen  somewhat  arbitrarily 
as  H  which  seems  to  be  a  suitable  interval  between  printouts  of  the  results. 


To  test  the  validity  of  the  model  we  can  compare  with  finite  dif¬ 
ference  calculations  in  one  dimension.  The  current  model  is  easily  modified 
to  planar  geometry  so  we  can  compare  for  both  spherical  and  planar  geometry. 


For  an  exponential  atmosphere  in  planar  geometry  we  have  the  results 
of  Enstrom  and  Brode7  as  well  as  a  local  Lagrangian  hydro-code  (HYDRO)  with 
which  to  compare.  In  Figure  11  two  cases  are  shown  using  these  two  methods. 

In  Figure  12  similar  results  are  shown  for  spherical  geometry. 

For  propagation  into  a  real  atmosphere  we  have  used  the  CIRA1  1  Model  5 
(hour  8)  case  giving  mid-solar  cycle  conditions.  Since  the  real  atmosphere 
is  far  from  isothermal  in  the  region  above  100  km,  one  expects  substantial 
differences.  In  Figure  13  planar  geometry  is  shown  for  the  current  model 
and  hydro-codes  while  Figure  14  gives  the  results  for  spherical  geometry. 

All  these  results  indicate  that  the  model  provides  a  reasonable  description 
of  the  propagation  so  long  as  the  relative  overpressure  is  much  less  than  100. 


The  model  is  readily  changed  for  propagation  at  some  angle  0 
relative  to  vertical  by  making  the  substitition  H  -*•  H  cos9  in  all  express¬ 
ions.  Having  done  this  we  can  compare  the  model  with  the  results  of 
Greene  and  Whitaker12  who  consider  the  initial  condition  of  a  4  MT  isothermal 
sphere  at  sea  level  using  both  a  one -dimensional  and  two-dimensional  hydro¬ 
code  for  propagation  at  45°.  The  model  was  given  initial  conditions  for 
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A-12.  Comparison  of  MODEL,  HYDRO  Code  and  Strong  Shock  Slope  fo 
Spherical  Exponential  Atmosphere. 
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Particle  Velocity  (km/sec) 


L0  and  So  using  Brode's  calculation  as  a  guide  as  discussed  in  the  following 
two  paragraphs.  The  results  are  shown  in  Figure  15  along  with  a  HYDRO  run 
with  the  same  initial  conditions  given  to  the  model .  Again  the  model  somewhat 
overestimates  the  pressures. 


The  model  makes  use  of  the  assumption  that  the  pulse  is  linear  in 
the  particle  speed.  To  test  this  assumption  a  HYDRO  calculation  was  made 
starting  with  a  linear  pulse  propagating  upwardly.  The  resulting  pulse  is 
shown  in  Figure  16  as  a  function  of  Lagrangian  and  Eulerian  coordinates. 

The  near  linearity  of  the  particle  speed  as  a  function  of  the  Lagrangian 
position  validates  the  assumption  of  linearity  for  the  model. 


We  are  now  prepared  to  use  this  model  for  a  variety  of  yields  and 
altitudes  to  determine  the  character  of  the  shocks  in  the  upper  atmosphere 
which  may  lead  to  communications  effects.  In  order  to  initialize  the  pulse 
we  have  used  Ihe  results  of  Brode  for  the  relative  overpressure  EC*)  and 
positive  pulse  duration  D*(r)  for  a  spherical  shock  wave  for  some  radii 
at  which  the  pulse  develops  a  linear  character.  For  radii  scaled  to  a 
dimensionless  variable  A  by 

‘/3 


r  -  x(^) 


EAR 


sc 


(A-69) 


we  take  the  overpressure  as 


g=  ^137  +  ^119  +  ^69 


.019 


(A-70) 


for  .  1  ?  6  ?  10  which  corresponds  to  2>A  Z  -25. 

Using  Brode's  D*(A)  and  taking  the  shock  speed  as  the  sound  speed 
for  this  weak  shock  region  the  pulse  length  was  fit  by 

L(A)  =  . 116A  +  .15  -  .0012  exp(3.84A).  (A-71) 


V 
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Height  (km) 

Figure  A-16.  Particle  Speed  in  a  Pulse  as  a  Function  of  Lagrangian  and  Eulerian 
Coordinates  as  Obtained  by  a  HYDRO  Calculation  for  an  Initially 
Linear  Pulse. 


In  order  to  initialize  0  and  L  for  all  angles  to  guarantee  equal 
energy  in  all  directions  the  initial  conditions  were  chosen  as 


Ro  =  Rsc 

L0  =  L(A)  x  R  (A- 72) 

sc 

and 

3  =  3(A) . 


H  T 

*  1.2  R  I 
sc  J 

This  procedure  sets  upper  limits  on  the  yield,  for  a  given  altitude,  which 
can  be  handled.  This  is 


Here  A  =  MIN  2 


and 


■*  dh) 

'  sc  / 


Jtn(x)  (cos0) 


Y 

max 


(A-73) 


which,  for  example,  means  Y  (sea  level)  =300  NTT  and  Y  (100  km)  =  .09  KT. 

max  max 

For  cases  outside  this  range  the  upward  traveling  shock  would  remain  very 
strong  at  all  altitudes  if  the  scale  height  were  constant.  In  Table  1  are 
shown  results  for  a  variety  of  yields  and  altitudes  for  propagation  upward 
and  for  0  =  45°. 


Sears  and  Kamm13  have  attempted  to  account  for  some  F-region 
disturbances  following  several  high  altitude  tests,  Teak  being  the  most 
prominent,  as  being  due  to  passage  of  the  shock  wave.  The  arrival  times  for 
the  disturbances  were  taken  from  ionosonde  data  and  shock  pressures  were 
calculated  based  on  these  times  and  empirical  scaling  laws.  They  concluded 
that  the  shock  waves  were  Mach  2  to  3  in  the  F-region  over  Maui  (1500  km  away 
from  the  Johnston  Island  test  site) . 


The  model  described  in  this  report  (MODEL)  would  not  generally  be 
applicable  to  the  cases  considered  by  Sears  and  Kamm.  However,  due  to  the 
nearly  horizontal  direction  of  propagation,  we  were  able  to  use  MODEL  with 
an  initial  shock  strength  consistent  with  the  required  pulse  shape.  It  was 
found  that  the  F-region  arrival  times  at  Maui  predicted  by  MODEL  were  2  to 
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Table  A-l .  Model  Results  for  Pulse  Relative  Overpressure  3  and 
Pulse  Length  L  for  Different  Yields  and  Altitudes. 


Altitude  (km) 

100 

150 

200 

250 

300 

350 

Yield 

(kT) 

Burst 

Altitude 

(km) 

Angle 

Degree] 

8 

L 

(km) 

6 

L 

(km) 

8 

L 

(km) 

B 

L 

(km) 

B 

L 

(km) 

B 

L 

(it  n) 

30000 

0 

0 

15. 

37 

13. 

115 

9.5 

150 

8.2 

180 

7.6 

210 

45 

7.5 

40 

8.0 

140 

6.2 

180 

5.4 

230 

5.2 

260 

10000 

0 

7.8 

27 

13. 

67 

8.1 

100 

6.3 

130 

5.5 

160 

5.3 

190 

45 

4.0 

28 

7r0 

70 

5.0 

115 

4.1 

150 

3.8 

190 

3.7 

220 

4000 

0 

4.5 

21 

8.0 

55 

5.5 

87 

4.5 

110 

4.1 

140 

4.0 

165 

45 

2.3 

21 

4.3 

58 

3.4 

95 

3.0 

130 

2.8 

160 

2.8 

195 

1000 

0 

2.1 

15 

3.8 

41 

3.1 

64 

2.8 

88 

2.6 

110 

2.6 

130 

45 

1.2 

14 

2.0 

36 

1.8 

69 

1.7 

95 

1.7 

120 

1.8 

145 

100 

0 

0.7 

7 

1.0 

20 

1.0 

33 

1.0 

45 

1.0 

60 

1.0 

70 

45 

0.4 

6 

0.6 

18 

0.6 

33 

0.6 

46 

0.6 

60 

0.6 

70 

2000 

33 

0 

5.0 

24 

7.4 

56 

5.0 

87 

4.0 

115 

3.6 

140 

3.6 

165 

45 

2.9 

28 

4.8 

67 

3.5 

105 

3.0 

140 

2.8 

170 

2.8 

200 

500 

0 

3.8 

22 

6.2 

52 

4.3 

80 

3.5 

105 

3.3 

130 

3.2 

155 

45 

2.0 

23 

3.3 

57 

2.6 

89 

2.3 

120 

2.2 

150 

2.2 

180 

100 

0 

1.6 

14 

2.6 

35 

2.1 

57 

1.9 

77 

1.9 

96 

1.9 

115 

45 

0  .9 

13 

1.4 

36 

1.2 

58 

1.2 

80 

1.2 

100 

1.2 

120 

10 

0 

0.6 

7 

0  .8 

18 

0.7 

30 

0.7 

40 

0.7 

50 

0.8 

60 

45 

0.3 

6 

0.5 

18 

0.4 

28 

0.4 

38 

0  .4 

50 

0-5 

60 

20 

67 

0 

1.7 

14 

1.8 

34 

1.4 

51 

1.3 

67 

1.3 

82 

1.3 

98 

45 

1.1 

16 

1.2 

38 

1.0 

57 

0  .9 

76 

0.9 

95 

1.0 

115 

5 

0 

1.5 

13 

1.5 

30 

1.2 

46 

1.1 

61 

1.1 

75 

1.2 

90 

45 

0.8 

13 

0.9 

30 

0.7 

47 

0.7 

62 

0.7 

77 

0.7 

92 

1 

0 

0.6 

7 

0.7 

18 

0.6 

28 

0.6 

37 

0.6 

46 

0.6 

55 

45 

0.4 

8 

0.4 

18 

0.4 

28 

0.4 

37 

0.4 

45 

0.4 

54 

.05 

100 

0 

- 

- 

0.4 

13 

0-3 

19 

0.3 

24 

0.3 

29 

0.3 

34 

45 

- 

- 

0.2 

14 

0-2 

19 

0-2 

23 

0.2 

27 

0.2 

32 

I 


3  times  those  observed  and  the  shock  pressures  were  S  to  10  times  less  than 
those  given  by  Sears  and  Kasim.  The  disparity  is  so  great  as  to  lead  us  to 
suspect  that  straight  line  shock  propagation  cannot  account  for  the  disturbances. 
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APPENDIX  B 

(From  an  unpublished  report) 

THE  CALCULATION  OF  LINEAR  DISTURBANCES  IN  AN  ISOTHERMAL  ATMOSPHERE 


The  simplest  model  of  the  atmosphere  which  exhibits  a  substantial 
number  of  the  features  of  interest  to  the  considerations  of  this  report  is 
that  of  an  isothermic  atmosphere  having  a  constant  ratio  of  specific  heats,  y, 
in  which  we  assume  all  disturbances  are  sufficiently  weak  that  the  linearized 
hydrodynamic  equations  apply.  If  we  restrict  ourselves  to  point  sources 
(which  depend  on  time) ,  we  can  obtain  an  exact  solution  with  no  great  effort 
or  expense.  We  have  written  a  small  computer  code  which  performs  this  task. 
The  reasons  for  writing  such  a  code  are  two-fold:  first,  the  model  in  hand 
approximates  the  earth's  atmosphere  and  the  pressure  source  due  to  an  under¬ 
ground  explosion  sufficiently  well  to  allow  useful  comparisons  with  test  data 
or  order-of-magnitude  estimates  of  some  parameters  in  hypothetical  situations; 
second,  we  can  use  the  exact  results  from  the  code  to  test  possible  approxi- 
maton  procedures  which  appear  potentially  useful  for  calculations  with  more 
complicated  atmospheres  or  sources. 

It  is  shown  in  Reference  B-l  that  the  response  of  an  isothermal 
atmosphere  to  a  point  source  at  the  origin  can  be  calculated  using  the  Green’s 
function 


G(z ,r,w) 


exp 


,  2  2.1/2 

(w  -wj 


/  2  2072 

|u>  -w  ) 

\  8/ 


(B-l) 


B-l 


.  -  AX.* 


where  z  is  the  altitude,  r  is  the  horizontal  distance,  and  tu  is  the 
Fourier  transform  variable  conjugate  to  the  time.  The  constants  in  Equation 
B-l  are  defined  as  follows 


a,  =  (y-l) 1/2  f 
g  c 


JA  2c  ’ 


,  2  2,1/2 

R  =  (r  +z  ) 


20) 

to  =  "'jj®'  and 
c  K 


R 

tO  =  c 


(B-2) 


where  g  is  the  constant  acceleration  of  gravity  and  c  the  speed  of  sound. 
We  can  use  this  Green's  function  to  calculate  the  relative  pressure  flucuation 
6p/po»  Pc  being  the  unperturbed  pressure,  by  computing 


Po 


(z,r,t) 


e2H  J  G(z,r,oj)  M(to)  e'1U)tdo) 
-00 


(B-3) 


where  M  is  the  Fourier  transform  of  the  pressure  source: 

00 

M(w)  =  J  M(t)eiWtdt  f 
-00 


(B-4) 


and  M(t)  is  the  relative  pressure  source.  As  it  stands,  Equation  B-3  is 
only  a  formal  expression,  for  we  have  not  chosen  the  branch  on  which  G  is  to 
be  evaluated  or  indicated  the  integration  path  (it  is  not,  in  general,  the 
real  axis).  From  (B-l),  we  see  that  G  has  branch  points  at  ±u»  ,  ±u>.,  and 

O  " 

±U)g.  These  points  are  real  and  their  relative  magnitude  is  as  shown  in 
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Figure  B-l.  We  choose  the  function  to  have  the  set  of  branch  cuts  shown  in 

Figure  B-l,  cuts  joining  ui.  to  ui  ,  ui  to  u  ,  hi  to  -«u  ,  -u>  to  -ui  ,  and  -uj 

A  g’  g  c»  c  c’  c  g’  g 

to  -  u»^,  all  the  cuts  running  along  the  real  axis.  For  real  values  of  u» 
greater  than  U)^,  we  choose  all  the  square-roots  in  Equation  B-l  to  be  evalu¬ 
ated  on  the  principal  branch  of  the  square-root  function.  With  the  chosen  set 
of  branch  cuts,  the  function  G  is  now  determined  everywhere  by  analystic 
continuation.  The  path  of  integration  for  Equation  B-3  is  that  shown  in 
Figure  B-l  (a)  running  slightly  above  the  real  axis.  That  this  path  must  be 
chosen  is  dictated  by  causality. 


For  most  choices  of  M(t) ,  the  integral  of  Equation  B-3  must  be  done 
numerically.  The  best  way  to  execute  such  a  calculation  will  depend  upon  the 
singularities  of  tf(ui)  and  the  behavior  of  this  function  for  large  values  of 
ui  .  Rather  than  attempt  a  general  discussion,  we  consider  the  class  of 
functions : 


M(t)  = 


(B-5) 


This  set  of  functions,  sometimes  called  N-waves,  provides  a  reasonable  repre¬ 
sentation  of  the  pressure  source  resulting  from  an  underground  explosion. 
Using  Equation  B-4,  we  calculate 


M(w)  = 


2Ai 

r  t 

+  1 

T 

-iute^ 
l  +  e  2 

1 

1  \ 

UI 

e 

2  iT  ui 

0 

'  C  ( 

2 

iTo 
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Figure  B-1.  (a)  The  location  of  the  branch  cuts  and 

integration  contour  defining  Equation  B-3; 
(b)  The  integration  contour  used  to  compute 
I-|  and  Ig. 


i 
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We  use  Equation  B-6  to  rewrite  Equation  B-3  as 


,  (B-7, 

*0 


where 


T 


_co 


(B-8) 


and 


-OO 


Looking  at  Equation  B-l,  we  see  that,  if  u>  has  a  large  imaginary 
part,  the  behavior  of  the  Green's  function  is  dominated  by  elu)to  .  The  be¬ 
havior  of  the  integrand  of  Equation  B-8  is  thus  dominated  by 
exp  iu)(tQ  +  Tq/2  -  t) ,  while  that  of  Equation  B-9  is  dominated  by 

exp  iu)(t  -  T  /2  -  t).  For  times  less  than  t  -  T  /2,  we  can  close  the  con- 
r  o  o  oo 

tours  for  both  1^  and  with  a  large  semi-circle  in  the  upper  half¬ 

plane.  Since  neither  G  nor  M  have  singularities  in  the  upper  half-plane, 
we  thus  get  zero  for  6p/pQ,  a  result  much  to  be  expected  since  tQ  -  To/2  is 
the  arrival  time  of  a  pulse  traveling  at  velocity  c.  If  t  lies  between 
tQ  -  Tq/2  and  tQ  +  Tq/2,  the  contour  for  1^  can  still  be  closed  in  the 
upper  half-plane  and  the  integral  again  gives  zero.  For  values  of  the  time  in 
this  range,  if  we  wish  to  close  the  contour  for  Ij  we  must  do  so  in  the  lower 
half-plane  and  thus  encircle  the  singularities  of  G.  ^  will,  therefore  be 
nonzero  for  tQ  -  Tq/2  <t.  In  the  same  way,  for  tQ  +  Tq/2  <t,  1^  will  be 
nonzero  and  both  1^  and  Ij  will  contribute. 
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As  discussed  in  the  previous  paragraph,  whenever  1^  or  is 
nonzero,  we  can  close  the  contour  in  the  lower  half -plane.  As  there  are  no 
singularities  in  the  lower  half-plane  or  on  the  real  asix  for  u»  >  uj^,  we  can 
contract  the  contour  to  be  the  one  shown  in  Figure  B-l(b).  It  may  appear  that 
we  can  contract  the  contour  further  to  integrate  the  discontinuity  across  the 
branch  cuts  along  the  real  axis.  This  procedure  would,  however,  lead  to  a 
problem  equivalent  to  integrating  exp(l/x)  through  the  origin  of  x.  We  must, 
in  fact,  choose  the  contour  with  some  care.  If  we  come  too  close  to  the  real 
axis,  the  singularities  there  cause  problems,  while  if  we  go  too  far  above  the 
real  axis,  the  exponential  increase  in  the  integrand  causes  difficulties. 
Using  a  simple  first  order  integration  technique  and  choosing  the  contour 
somewhat  carefully,  we  can  achieve  three  figure  accuracy  by  using  a  few  thou¬ 
sand  points. 


If  u).  =  u>  (corresonding  to  a  value  y  =  2),  the  worst  singulari- 
ties  in  the  Green's  function  cancel,  and  we  can  indeed  contract  the  contour 
all  the  way  to  the  real  axis  since  the  discontinuity  of  the  integrand  in  that 
case  contains  only  integrable  square-root  type  singularities.  We  can  thus 


save  the  trouble  of  choosing  a  contour,  and  the  approximation 


u>.  =  U) 

A  8 


is  one 


whose  accuracy  we  wish  to  investigate.  Another  simplifying  approximation, 


called  the  low  frequency  limit,  consists  of  assuming  that  the  only  important 


contributions  to  the  integrals  Ij  and  1^  come  from  the  branch  cut  from  u»c 
to  -u>£ .  In  this  approximation  an  answer  can  be  given  in  closed  form.  In 
Reference  B-l,  several  examples  are  given  in  which  the  low  frequency  limit 


gives  very  accurate  answers  from  very  soon  after  t  on.  It  appears  that  this 
good  agreement  is  not  general,  however,  but  results  from  the  particular  way 


the  source  term  used  in  Reference  B-l  depends  on  the  time. 


As  an  illustration  of  the  results  of  the  calculations,  we  consider 
the  particular  set  of  parameters: 
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uig  =  1.049  x  10  2  sec  1 

u).  =  1.068  x  io  2  sec  * 
A 

T  =  100  sec 


o 

z 

r 

c 


2  x  10”  cm 
n 

,4 


107  cm 


=  3  x  io  cm/sec 


We  choose  the  amplitude  of  the  initial  disturbance  such  that  the 

amplitude  of  the  N-wave  arriving  at  the  point  (z,r)  will  be  1.  In  Figure  B-2 

we  show  the  results  for  times  later  than  450  seconds  after  the  source  first 

became  nonzero  (the  starting  time  for  the  plot  was  arbitrarily  chosen;  the 

result  is,  of  course,  nonzero  from  t  -T  /2  on).  Plotted  on  the  large  graph 

oo 

are  the  results  for  setting  u»^  and  u»^  both  equal  to  1.049  x  10  and  for 
the  low  frequency  limit.  The  low  frequency  limit  completely  misses  the  early 
time  structure  of  the  result.  This  structure,  which  is  substantial,  might  be 
important  to  considerations  of  detection  or  coupling  to  other  atmospheric 
modes.  The  result  of  setting  u)^  =  appears  quite  satisfactory  in  this 
particular  case,  being  very  close  to  the  exact  result  at  all  times. 


In  Figure  B-3,  we  show  the  results  of  the  case  of: 


u>  =  8.25  x  io'3  sec"1 
8  -3  -1 

u).  =  8.75  x  10  J  sec 


=  19.4  sec 

=  2.25  x  io7  cm 
=  1.3  x  107  cm 

A 

=  4.3  x  10  cm/sec 


We  show  the  results  for  times  later  than  900  seconds;  the  earlier 

results  are  off  the  graph.  The  results  are  shown  for  the  exact  calculation 

-3  -1 

and  for  u».  =  u>  =  8.25  x  io  sec  .  The  results  of  the  low  frequency  limit 

A  g 
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